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ABSTRACT 

We present a numerical study of the global modes of oscillation of thick accretion discs 
around black holes. We have previously studied the case of constant distributions of 
specific angular momentum. In this second paper, we investigate (i) how the size 
of the disc affects the oscillation eigenfrequencies, and (ii) the effect of power-law 
distributions of angular momentum on the oscillations. In particular, we compare 
the oscillations of the disc with the epicyclic eigenfrequencies of a test particle with 
different angular momentum distributions orbiting around the central object. We find 
that there is a frequency shift away from the epicyclic eigenfrequency of the test 
particle to lower values as the size of the tori is increased. We have also studied 
the response of a thick accretion disc to a localized external perturbation using non 
constant specific angular momentum distributions within the disc. We find that in 
this case it is also possible (as reported previously for constant angular momentum 
distributions) to efficiently excite internal modes of oscillation. In fact we show here 
that the local perturbations excite global oscillations (acoustic p modes) closely related 
to the epicyclic oscillations of test particles. Our results are particularly relevant in 
the context of low mass X-ray binaries and microquasars, and the high frequency 
Quasi-Periodic Oscillations (QPOs) observed in them. Our computations make use 
of a Smooth Particle Hydrodynamics (SPH) code in azimuthal symmetry, and use a 
gravitational potential that mimics the effects of strong gravity. 

Key words: accretion discs — black hole physics — hydrodynamics — stars: neutron 
— X-rays: binaries 



1 INTRODUCTION 

Thick accretion discs may appear in systems containing 
compact objects such as black holes and neutron stars. These 
could form in different scenarios governed by gravitational 
forces. Specific examples are: (i) gravitational collapse of a 
rotating massive star, (ii) the coalescence of two compact ob- 
jects or (iii) accretion phenomena asociated with Low Mass 
X-ray Binaries (LMXBs). 

The first case (gravitational collapse), may be related 
to the production of 7-ray bursts (GRBs) as in the collap- 
sar or hypernova models. The co llapsar (or failed supernova) 
was proposed bv lWooslevI il993Tl and offers the possibility to 
form a thick disc when the mantle of a rotating massive star 
(M ~ 25Mq ) falls freely onto the recently formed black hole. 
The material has a very high angular momentum, which pre- 
vents it from directly accreting and a centrifugally supported 
disc is formed. These systems are related to energetic, long 
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durat ion (> 2 s) GRBs accompanied by a Type Ib/Ic super- 
nova (iKulkarni et alJll998t IStanek et aT]l2003l : iHiorth et all 

2003). Similar numerical relativistic calc ulations concerning 
the co llapse of massive stars, reported bv lShapiro fc Shibatal 
(2002) and Shi bata &i Shapiro! I|2002j) may explain the for- 
mation of massive black holes surounded by a debris torus . 
The hypernova scenario was proposed bv iPaczvriskil jl998t) 
and is also related to GRBs. Hypernova are extremely 
energetic supernova explosions (kinetic energy ~ 3 x 10 52 
erg), which consider (a) a single star (i.e. Wolf-Rayet 
type star) or (b) a binary system formed by massive stars 
(M ~ 20 - 25Af ) orbiting closely (~ 1.5Rq) around each 
other, and which eventually merge to form a single star [see 
iNomoto. Maeda. Mazzali. Umeda. Deng fe Iwamotol J2004f) 
for a comprehensive review] . This star will undergo gravita- 
tional collapse, leading to the formation of a compact object. 
Due to the rapid rotation, a thick accretion disc around the 
remaining black hole will be formed. In either hot, 
thick disc of ~ O.IM© is expected around the newborn black 
hole. 
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Another posibility to form thick discs is via the coales- 
cence of two compact objects, leading to the formation of 
a system composed of a black hole and an accretion torus. 
This has been investigated by various groups using num- 
m erical simulations. Some of the first r esults were obtained 
bv lDavies. Benz. Piran fc Thielemann! l)l994l) . who modeled 
the merger of two neutron stars. They showed that after few 
orbital periods the residual material could form a thick disc. 
These results have been confirmed and extended (see e.g., 
iRuffert. Janka fc Schaefe3ll996t iRuffert fc Jankalll999l) for 
the merger of two neutron stars, where again, a hot thick disc 
is formed around a central mass, which will presumably col- 
lapse to a black hole. Investigating t he more complex fusion 
of a neutron st ar and a black hole, iKluzniak fc Leel ill 9981) 
and iLeei {2001) found that after the initial mass transfer 
there is a remnant of material around the black hole (see 
also recent results bv lRosswog. Speith fc Wvnnll20(5) . All 
these scenarios may form massive (Md/MsH ~ 0.1) and 
dense discs. 

There are also low-mass discs (Ma/ Mbh <C 1) which 
are formed through mass transfer in a binary system. A 
steady supply of mass, energy and angular momentum 
produces structures which are nearly in a steady state, but 
which ne vertheless e xhibit variability on various timescales 
(see e.g.. lKatolll99l) . 

Once the disc is formed, one may ask: under which con- 
ditions is it dynamically stable or unstable (quite apart from 
the question of thermal stability, if dynamically stable)? An- 
swering this in any detail requires the consideration of vis- 
cosity, self-gravity, magnetic fields and the angular momen- 
tum dist ribution of the material wit hin the disc. This last 
point led lPapaloizou fc PringlJ (Il984f) to investigate the sta- 
bility of an isentropic and constant angular momentum disc, 
and find that it is unstable to global non-axis ymetric pertur- 
bation s. Three-dimensional calculations by IZurek fc Bend 
(1986) found that there is indeed a redistribution of angu- 
lar momentum within the disc, due to the growth of non- 
axisymetric instabilities. In the presence of accretion, how- 
ever, this instability may be suppressed feraeslll987l) . and 
thick tori with relatively flat distributions of angular mo- 
mentum may actually occur in astrophysical systems. 

Perhaps one of the most interesting global instabilities 
is the so-cal led runaway radial instability , whic h was 
identified by lAbramowicz. Calvani fc Nobilil l|l983l) . It is 
due to a violent mass exchange asociated with the rapidly 
changing equipotential Roche surfaces of the disc + black 
hole system. It is catastrophic and leads to the complete 
destruction of the disc on a dynamical time scale. It is well 
known that in binary systems there is a point called the 
interior Lagrangian point L\, located at a radius rz, x where 
there is a balance between gravitational and centrifugal 
forces. For accreting discs, when mass transfer from the disc 
to the central object begins the corresponding Lagrange 
point, lying between the compact object and the disc, 
moves outward. The accretion r ate incr e ases r apidly until 
the disc is accreted completely. IWilsonl il984T) found that 
the rotation of the bl ack hole inhibits the runaway radia l 
instability and later lAbramowicz. Karas fc Lanzal Jl998h 
confirmed this result and reported that the self-gravity 
of the disc has a destabilizing effect. It appears that the 
most critical parameter in the runaway radial instability 



is the specific angular momentum distribution within the 
disc. In the last twenty five years there has been much 
work on this issue from stationary (i.e. Schwarzschild 
metric) to dynamical (i.e. Kerr metric) space-time models. 
A complete su mmary of the work d one appears in the 
extens ive work of lFont fc Daignd <2002ft andlDaigne fc Font] 
(2004). These authors conclude, as iDaignTfcMo^hkovitcI] 
(1997) did earlier, that the instability is suppresed when 
the angular momentum follows a power law distribution 
with a positive radial gradien t. I n this resp ect studies 
made bv hluffert fc JankaNl999l) andlEeel (120011) concerning 
compact binary mergers show that the distribution of 
angular momentum in post-merger discs is far from being 
constant, and thus they are stable in this context. 

Our numerical simulations were performed in order to 
address two questions: (i) How do thick discs oscillate? and 
(ii) How does a localized perturbation affect the global os- 
cillatory behaviour of the disc? These are formulated in 
the particular context in the millisecond oscillations discov - 
ered in X-ray binary sources by RXTE (Ivan der Klislfc oOO). 
These systems have compact objects surrounded by an ac- 
cretion disc, whose variability could be imprinted on the 
X-ray lightcurve and corresponding power spectra. 

Disc oscillations have been investigated numerically 
and analytically by several groups before, and may be ap- 
plied to the quasi perio dical oscillations (QPOs) observed 
in many X-ray sourc es ||Perez_et_alJ 11997J; |Waronei| Il999l : 
ISilbergleit et al.lll999l: lOrtega-Rodriguez et al.ll2002l). One 
specifi c model llRezzc41^^^shida^Maroarone" fc ZanottH 
2003a) accounts for the high frequency QPOs as a result 
of p-mode oscillat ions of an accretion torus orbiting c lose 
to the black hole. iRezzolla. Yoshida fc Zanottil (l2003bl) in- 
vestigated theoretically how a relativistic torus oscillates 
in a stationary Schwarzschild space-time, showing that the 
eigenfrequencies of the axisymmetric oscillations correspond 
to acoustic p-modes. Their results shows that the oscilla- 
tions behaves like a sound wave globally traped within the 
disc, with the eigenfrequencies appearing in the sequence 
2:3:4... independently of the distribution of angular mo- 
mentum considered. The same con c lusion s are obtained by 
IZanotti. Font. Rezzolla fc Monterol ll2005l) using numerical 
simulations for a dynamical Kerr space-time, concluding 
that p-mode oscillations in relativistic tori could explain the 
high frequency QPOs observed in the X-ray binaries. The 
oscillatory analysis made by Rezzolla et al. (2003a) , Rezzolla 
et al. (2003b), Montero et al. (2004) and Zanotti et al. (2005) 
has shown theoretically that there is a frequency shifting to 
lower frequencies as the size of the torus is increased. Their 
analytical work is restricted to height -integrated structures 
and globally applied perturbations. We have confirmed this 
result more generally here with a different approach. 

Fro m the observations in th e X-ra y band mentioned 
above, lAbramowicz fc Kluzniakl i200lf) noted in GRO 
J1655-40, that the two stable peaks that appear in the power 
spectrum in the hHz range are in a 3:2 ratio. The same com- 
measurability has now been observed in three more sources: 
H1743-322, XTE J1550-564 and GRS 1915+105 (see 
iMcClintock fc Remillardl2004 for a recent review). This be- 
haviour may also be present in a system containing a neu- 
tron star (Sco X-l) llAbramowicz. Bulik. Bursa fc Kluzniakl 
2003), and could be explained by parametric resonance 
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in a thin disc iAbramowicz et all I2QQ3I : iRebuscol |2004) . 
These models predict that the frequencies will show the 
epicyclic motions of a perturbed flow line in the accre- 
tion disc or a combination of these and a fixed pertur- 
bation associated to the spin of the central object (in 
the neutron star sources). This last issue was analysed by 
iKluzniak et all ll2004).lLee Abramowicz fc Kluzniakl (120041) 
and lRublo^Herrer^fcLe^T2005l hereafter Paper I) showing 
that it is indeed posible to excite internal oscillatory modes 
efficiently by an external driving agent. 

In this paper we show that there is a shifting to lower 
frequencies as the thickness of the disc is increased, inde- 
pendently of the distribution of angular momentum within 
the torus. We also show that it is possible to excite internal 
p-modes in a thick accretion disc, by means of a localized 
perturbation which affects only a small portion of the disc 
periodically, exciting the strongest modes apparent in a 3:2 
frequency ratio. 

We empasize that the results presented here are purely 
dynamical in a pseudo-Newtonian gravitational potential 
which mimics the properties of the Schwarzschild metric 
for a stationary space-time. No attempt has been made to 
consider magnetic fields, radiation or the effects of viscos- 
ity. For definiteness, a numerical value of 2.5Mq has been 
used throughout the paper. The reader should keep in mind, 
though, that distances are always measured in units of the 
gravitational radius. As such, our results are scale free, and 
can be applied to any black hole mass, with the frequencies 
scaling as v oc 1/M. 



2 INITIAL CONDITIONS AND NUMERICAL 
METHOD 

2.1 Hydrostatic equilibrium for a thick torus 

We will define a thick torus as that which has similar radial 
and vertical extensions, i.e., L ~ H, and where L is compa- 
rable to its radial position, R. We will neglect self-gravity 
and additionally assume azimuthal symmetry. We have 
constructed thick tori using th e pseudo-Newtonian grav i- 
tational potential proposed bv IPaczvriski fc Wiital |l980), 
given by: 

-GMbh 



$ = 



R 



(1) 



which qualitatively reproduces the behaviour of a test par- 
ticle in the Schwarzschild metric, and accounts correctly for 
the positions of the marginally stable and marginally bound 
orbits. In this equation, R — \/r 2 + z 2 is the radial distance 
mesured from the central mass and r g — IGMbh /c 2 de- 
notes the Schwarzschild radius of the central object. Using 
a polytropic equation of state P = Kp 1 , the equations of hy- 
drodynamics in cylindrical coordinates, and neglecting self 
gravity, we have: 

1. 



-VP: 



-V$ c 



(2) 



where <J> e g is the effective potential. Using the Paczyriski- 
Wiita expression given in (1), substituting it in (2) and in- 
tegrating over r we have, 



7 



7-1 P 



+ $ c ff + $ = est, 



(3) 





Figure 1. Meridional cross-sections over one half the i — z plane 
for tori with constant specific angular momentum, I = const 
(above), and a power-law specific angular momentum distribu- 
tion I oc r a , a = 0.1 (bottom). The filling factor $o is indicated 
in units of c 2 /2. 



which describes a torus in hydrostatic equilibrium. The term 
<E>o can be interpreted as a filling factor of the effective po- 
tential well, and will be discussed below. Using an arbitrary 
distribution of angular momentum £(r), which depends only 
of the position coordinate r according to the von Zeipel theo- 
rem for a polytropic equation of state, the effective potential 
takes the form: 



-GMbi 
R-r a 



+ 



2r' 3 



■dr. 



(4) 



The boundary of the tours is the surface over which the 
pressure is zero, and is given by: 

1/2 

J\2 



GMbi 



£(r' 



2r' 3 



+ 3>o 



(5) 



The corresponding meridional cross sections over one half 
the 1 — z plane are shown in Figure 1 for various filling fac- 
tors $o and different distributions of specific angular mo- 
mentum. Several groups have studied the dynamics and 
stability of tori with varying distributions of specific angu- 
lar momentum jDaigne fc Fondl2004t iFont fc Daignell2002l: 
[ Masuda. Nishida fc Eriguchil Il998t iDaigne fc Mochkovitcbl 
1997). In summary, the conclusion has been that in a torus 
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Figure 2. Angular momentum distributions used for our models. 
The solid line represents a Keplerian angular momentum distri- 
bution for test particles in circular orbits, while the dashed line 
shows the non-constant power law distribution. For the case with 
uniform angular momentum we used i = 1.92 r g c. 



with non-constant angular momentum, the runaway radial 
instability is suppresssed. We have performed calculations 
(described below) in this context, and have used the value 
£ = 1.92r g c for a constant distribution. This corresponds to 
the angular momentum of a test particle in a Keplerian mo- 
tion at 
used 



r = 4.25r 9 . For the power law distribution case we 



where 



1/2 


'GM BH r 3 ' 


/ r=r 


(r - rg) 2 _ 



1/2 



(6) 



(7) 



The value of the power index a, can be super-keplerian 
a > 0.5, keplerian a = 0.5 or sub -keplerian a < 0.5. T o 
compare with the results obtained bv lFont fc Daignd (2002), 
we have used a sub-keplerian distribution with an index 
a = 0.1. Finally the normalizing factor for the radius is that 
which corresponds to the Keplerian value ro = 5.1r s . We 
will henceforth define the centre of the torus as the locus (a 
circle, really) of maximum density and pressure in the disc. 



2.2 Numerical method 

We have used the SPH (Smooth Particle Hydr odynamics) 
nume rical method to perform our simulations jMonaghanl 
1992). Essentially this method works by interpolating a de- 
sired function A(r), with a set of given quantities A(r') 
throught a kernel function W(r, h), using the following con- 
volution integral: 



A(r) = / A{r')W{r-r')dr. 



(8) 



The implementation of t he code was made usi n g the 
prescriptions calculated by iMonaehan fc Lattanzid (Il985h 
when the system has azimuthal symmetry, using the follow- 
ing kernel function: 



W(r, h) 



a 

1^ 



l-3(r/ft) 2 /2- 

{2-r/hf/A 





■3(r/ft) 3 /4 



if < r/h < 1; 
if 1 sC r/h < 2; 
2 sC r/h. 



In the last equation, h represents the smoothing length, and 
is comparable to the typical separation between fluid ele- 
ments, and r is the radial coordinate (h is essentially the 
spatial resolution of the calculation). As our calculations 
were performed in two dimensions, v — 2 and a = 10/77;. 
We did not include the viscosity of the gas and hence the 
Navier-Stokes equations of motion take the form: 

dv r IdP GM B Hr , ^ 2 , (dv r \ 

\ dt J art 



dt 
dv- 
~~dt 



R(R-r g y 

GM BH z 
R(R-r g y 



(10) 



+ 



\ at J ar t 



(11) 



p dr 

1 dP 

p dz 

where again R — \Jr 2 + z 2 is the distance to the central 
object. The sub-index art indicates the artificial viscosity 
terms, used in the numerical calculations to account for the 
presence of shocks and avoid particle interpenetration. We 
have additionally for the conservation of angular momen- 
tum: 

dt 

And finally, for the energy equation we have: 



(12) 



de 
dt 



P 



v + 



V dt) ar t 



(13) 



where e is the internal energy per unit mass. We have as- 
sumed that we do not have losses of energy by external 
causes (i.e. AQ = 0) and the only changes in e arise from 
mechanical work done on the system. In discretized form 
for numerical work, the convolution integral given in equa- 
tion (8) becomes a sum over all the fluid elements in the 
calculation. 

For the artificial vi scosity, we have used the prescription 
given bv lBalsaral l)l995h which reduces the artificial shearing 
stress during the time evolution of the disc. The explicit 
forms for the acceleration and the energy dissipation due to 
artificial viscosity for the i — eth SPH particle are then: 



and 



mjTIyV iWi- 



V dt J i art 



2 ^ 



mjUij(vi — Vj)- ViWi- 



(14) 



(15) 



where, a gain, W is the kernel funct ion and n is defined by 
(see e.g.. lLee fc Ramirez-Ruizll2002ri 



= ( -4 + -§ I = {-ctbPij + p b p\ 



(16) 



( 



(gi-gjHfi-fj) (/,+/,) 

h i] (\r i -r- j \yhf.) + v l 2c zj 



< 0: 



if (vi — Vj) ■ (fi~ 

if (vi — Vj) ■ (fi — fj) > 0; 

In this last equation the function fi is defined by: 

f = \y-vU fl „v 

Jl |V-% + |V xv\i + ri'ci/hi' K ' 

The factor rj ~ 10 -4 in the denominator prevents numeri- 
cal divergences, Ci represents the sound speed and hi is the 
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lenght scale for the i — eth particle, at = = 7/2 are 
constants of order unity and 7 is the polytropic index from 
the equation of state. This form of the artificial viscosity 
supresses the shearing stresses when the compression in the 
fluid is low and the vorticity is high |V • v\ -C |V X v\, but 
remains in effect if the compression dominates in the flow 
IV • v\ > IV x v\. 



2.3 Initial conditions 

Our initial conditions are generated by fixing the mass of 
the central object, the equation of state, with 7 = 4/3 and 
the constant K, the filling factor for the potential, <3?o and 
the distribution of angular momentum £(r). We then gen- 
erate the fluid elements which will model the torus using a 
Monte Carlo random number generator. The particles that 
fall inside the torus density profile will be the fluid elements 
of the simulation. These fluid elements are then relaxed for 
several dynamical times with an artificial damping term in- 
cluded in the equations of motion in order to obtain the a 
distribution which is as close to equilibrium as possible. We 
are able to build tori of different sizes by simply varying the 
filling factor $0 as was shown in Figure 1. 

Our initial configurations are always non-accreting tori, 
confined inside their Roche lobe, i.e., r; n > rz, 1 where ri n 
is the distance from the central object to the inner edge of 
the torus and tl x is the radius of the inner Lagrange point, 
L\. When we focus on the pure oscillation dynamics of the 
disc, we also wish to avoid accretion. In those cases, this 
is accomplished by maintaining the mass of the black hole 
constant (even if some particles overflow the Roche lobe as 
a result of the applied perturbation). 

Our code produces different snapshots of the main hy- 
drodynamical variables of the system like the position of the 
centre of the disc, and also the maximum and mean densities 
and the potential, kinetic and internal energies as functions 
of time. This data then is analysed by calculating the corre- 
sponding Fourier transforms and looking for characteristic 
frequencies. 



2.4 Testing the code 

To test our code, we perform ed simulations reproducing pre- 
viousl y publishe d resutl s jAbramowicz . Calvan i fc Nobiul 
' 19831 ). iDaiene fc Mochkoyitchl Jl997t) . 



lAbramow icz. Ka ras fc Lanzal il998f) . lFont fc Daignel ll200L. 
and IDaiene fc Fontl ll2004h 'l investigating the behavior of 
the runaway-radial instability with changing distributions 
of angular momentum (in systems without self-gravity). 
We constructed tori that (i) almost fill their effective 
potential wells, i.e., Ti n ~ vl x , and (ii) have constant and 
non constant angular momentum distributions. We then 
introduced radial impulsive perturbations, adding a velocity 
field to the initial di stribution. Th is field is similar to the 
solution proposed bv lMichel (Il972f) for spherical relativistic 
accretion, and has the form 



-V- y/r g /R ; . 



(19) 
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Figure 3. Mass in the disc as a function of time for two test 
runs. The dashed line corresponds to the model with constant 
angular momentum distribution within the disc (t, = est), while 
the solid line corresponds to the model with a power law angular 
momentum distribution i oc r a . 



We have found results that are in complete agreement 
with the previous work available in the literature, that is, the 
runaway instability is completly supressed when one consid- 
ers a non-constant angular momentum distribution within 
the torus. This is shown in Figure 3, where the mass of the 
disc is plotted as a function of time for a calculation with a 
power law distribution of specific angular momentum (solid 
line), and another with constant specific angular momentum 
(dashed line). The former is clearly stable, while the latter 
is not, and shows the characteristic runaway behavior. Once 
the instability is triggered, the torus is destroyed in a few 
dynamical times. 



3 NATURAL OSCILLATIONS OF THE DISC 

To study the free oscillations of the disc, we have performed 
calculations where small, impulsive perturbations have been 
applied to the velocity field, as in equation (19), only in the 
initial conditions. This was done in particular to study how 
the thickness of the disc affects its modes of oscillation. The 
space-time around the black hole was kept constant during 
this set of runs, even if slight overflow occurred through the 
inner Lagrange point. 



3.1 Epicyclic frequencies 

The results obtained were compared with the behaviour of a 
test particle in a circular orbit around a central mass. If we 
introduce a small radial perturbation, the particle will show 
radial oscillations at the epicyclic frequency k. Our study 
begins with the equation for radial motion of the particle: 



(20) 



where the parameter 77 <C 1 modulates the strength of the 
perturbation. 



(fr _ _d$off 
dt 2 ~ dr ' 

Introducing a perturbation of the form: Ar = r — ro, and 
calculating a first order Taylor expansion over r for small 
perturbations, we obtain 
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Figure 4. Epicyclic frequency function of radius. The 

curve is computed for a test particle orbiting a central object of 
2.5Mq. Compare with Figures 5 and 6. 



d 2 Ar 
dt 2 



dr 2 



Ar, 



(21) 



for a fixed radius. This equation is simply that of a harmonic 
oscillator, with epicyclic frequency k = (d 2 $ e g /dr 2 )^ . For 
a constant angular momentum (l = est), we have: 



n(r) 



1 GM BH (r-3r g 



2tt 



r(r 



(22) 



In Figure 4 we show the radial dependency of k(t). 



3.2 Results 

We show here results for tori of three different sizes: small 



( r d < r g ), intermediate (rd 



and large (r d > 



both for I =cst and £(r) = r a . In each case we com- 
puted the Fourier transform of the radial motions of the 
centre of the disc. This always exhibited one prominent 
peak, at frequency k*. In the limit of small tori, k* — * 
/t(ro), the value for test particle motion. As the size of 
the torus increases, there is a shift to lower frequencies, 
and k* decreases. The Fourier transforms and the corre- 
sponding values for the frequencies (along with the model 
parameters) are given in Figures 5 and 6, and in Ta- 
ble _l_ ; _Tlwse_jiumen^ 

byjRezzolla :i _Yoslii^ ^ d2003al) and by 

iRezzolla. Yoshida fc Zanottil i2003bT) . where the frequency 
shift has been computed analytically in the Schwarschild 
metric (albeit with the aforementioned simplification of 
height-integrated equations). We thus confirm the p-mode 
interpretation of Rezzolla et al. (2003), in which the radial 
epicyclic frequency at the position of the centre of the torus 
represents the value at which the fundamental p-mode fre- 
quency tends to in the limit of vanishing torus size. 



4 LOCALIZED PERIODIC PERTURBATIONS 
AND GLOBAL RESPONSE 

4.1 A localized external perturbation 

Having identified the strongest response in the radial os- 
cillations as a modified epicyclic frequency, or p-mode, we 
proceeded to apply localized, periodic perturbations to thick 
tori (as in Paper I). We believe it is reasonable to assume 
that the perturbation will be stronger in the inner regions 
of the disc, if it arises from the central object, and periodic, 
if it is related to its spin period (clearly identifiable in neu- 
tron star systems) or some other mechanism repeating at 
intervals At = l/u ver . This perturbation can be related to a 
wide variety of phenomena coming from the central object 
or from the material that conforms the disc. A specific ex- 
ample could be the emission of gravitational radiation from 
the black hole when a clump of material is accreted and 
breaks the axial symmetry in a cuasi-periodic fashion. This 
could produc e a time variation in the mass quadrupole of 
the disc itself dZanotti. Rezzolla fc Fontl2003t). Another pos - 
sibility in this context was proposed by Ivan Putter] fcOOll) . 
and is related to the connection between the magnetic field 
of the disc and the gravitational radiation produced by the 
black hole. Disc oscillations can be related to the forma- 
tio n of the disc, or to its intrinsic properties . Simulations 
by llgumenshchev. Chen fc Abramowiczl il99rJl have shown 
that varia blity and os cillations can be triggerred by viscous 
processes. iKatol i200 lT) has suggested that energy flux within 
the disc through viscous and angular momentum processes 
can alter the equilibrium state, producing oscillations. Fi- 
nally, in highly dynamical situations, the disc will naturally 
lack complete azimuthal symmetry (e.g., in compact object 
mergers) and the formation of clumps inside the disc could 
lead to small oscillations. 

We have modeled perturbations in the context of non- 
constant angular momentum discs in the same way that we 
did in Paper I, introducing an additional acceleration in the 
equations of motion, given by: 



Gtpert = —r\a g ■ exp 



/ r-p - r \ 
V Sr J 



Sm(2TTV per t). 



(23) 



Here a g is the acceleration due to gravity, ro is the outer 
edge of the torus and )) < 1 is a parameter that mod- 
ulates the strength of the perturbation. The exponential 
term decays on a scale Sr ~ R, the radial extension of the 
disc, thus reproducing the desired behaviour for the per- 
turbative force, which will be strong near the inner radius 
and weak in the outer regions. This acceleration induces ra- 
dial oscillations, which can be Fourier analyzed to extract 
the main frequencies in a similar f ashion as was done by 
iLee. Abramowicz fc Kluzniakl i2004l) . 



4.2 Results 

We performed simulations using constant and non-constant 
angular momentum distributions within the disc. For the 
constant distribution, we find, as described in Paper I for a 
few particular cases, that the localized perturbation induces 
global oscillations at various characteristic frequencies. The 
most prominent (besides that ocurring at the perturbing fre- 
quency Vper) is at V\ = k*, the modified epicyclic frequency 
associated with the radial motions of the centre of the disc. 
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Table 1. Models for discs of different size orbiting a 2.5Mq black hole. We give the extent of the 
disc r d , the disc to black hole mass ratio (Mj/Mbh), the angular momentum value (£o/r g c) or 
distribution (I oc r"), the locus of the point of maximum density (R c /r g ), the filling factor which 
gives us the size of the torus ($o/(c 2 /2)), the eigenfrcquency re* at which the point of maximum 
density oscillates, and the number of SPH particles we used in the simulation (TV). 



Model 


size (r d ) 


M d /M BH 


£o/r g c 


R c /r g 


*o/(c 2 /2) 


k* (Hz) 


N 


(a) 


r d < r g 


3.0 x icr 8 


1.9197 


4.25 


-0.1183 


425 


2337 


(b) 


r d ~ r g 


5.9 x 10~ 3 


1.9197 


4.25 


-0.1060 


355 


2886 


(c) 


r d > r g 


0.14 


1.9197 


4.25 


-0.0841 


290 


4082 


(d) 


r d < r g 


1.6 x 10~ 8 


oc r" 


5.10 


-0.0749 


350 


2317 


(e) 


r d ~ r g 


2.0 x 10~ 3 


ocr" 


5.10 


-0.0678 


300 


1866 


(f) 


r d > r g 


0.3 


oc r" 


5.10 


-0.0595 


250 


1630 



model (a) 




model (e) 





Figure 5. Meridional cross sections for a small, intermediate and large torus (left panels), for a 
constant angular momentum distribution (£ = est), and Fourier transforms of the radial oscilla- 
tions of the centre of the torus for each case (hereafter represented by the letters CMR in the 
figures) in the right panncls. The discs orbit a 2.5Mq black hole. Note how the eigenfrcquency 
re* is progressively shifted to lower values as the size of the torus increases and compare with the 
epicyclic eigenfrcquency for a test particle shown in Figure 4. For all the unperturbed tori shown 
here, the centre is located at r = 4.25r 9 . 
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model (d) 




7 4 6 8 10 

r/r, 



Figure 6. Same as Figure 5, but for a power— law 
all the unperturbed tori shown here, the centre i 

A second broad peak with about ten times less power at 
i/2 is in a 3:2 relation with the former (see Figure 7). This 
latter feature is barely made out in the calculations where 
an impulsive perturbation was applied to the whole disc at 
t = 0. 

For a power law distribution of angular momentum we 
find that the disc responds somewhat differently to the lo- 
calized perturbation. A typical power spectrum for the ra- 
dial motions of the centre of the disc is shown in Figure 8 
(when v per = 100 Hz). Two strong narrow peaks are evi- 
dent at v p£r and 2u ver , despite the fact that the perturba- 
tion is a pure sinusoid. Furthermore, a weak subharmonic 
at v per /2 is also present. This is clear evidence of non- 
linear coupling between various modes, and is reminiscent of 
analogous results for forced oscillations in slender tori pre- 
viously reported in the context o f kHz QPOs in LMXBs 
llLee. Abramowicz fc Kluzniakll2004) . Note that the over- 
tone is absent in the power spectrum in Figure 7 (when 
the angular momentum is constant in the disc). As for the 
case with constant angular momentum, two more peaks, at 
V\ = k* = 250 Hz, and V2 ~ 1.5vt can be seen in the 




?G0 400 600 BOO 1000 

Frequency (Hz) 



angular momentum distribution (£ oc r a ). For 
located at r = 5.1r 9 . 

spectrum. To the limit of our resolution, they are in a 3:2 
correspondence. Thus it appears that the angular momen- 
tum distribution is somehow responsible for a certain mode 
coupling which allows higher harmonics and subharmonics 
of the perturbation to appear in the power spectrum. Note 
also that when an impulsive perturbation was used in a thick 
torus with a power law distribution of angular momentum, 
secondary features, consistent with being in a 3:2 relation 
with the stronger ones, are apparent in the spectrum once 
the torus becomes relatively thick(see Figure 6, middle and 
bottom right panels). Apparently the periodic nature of the 
forcing also has a direct influence on the magnitude of the 
response of the system. We elaborate further on this point 
below in § 



5 DISCUSSION 

We have performed a numerical study of the global oscil- 
lations of geometrically thick fluid tori in axisymmetry. We 
have neglected self-gravity, but used a pseudo-potential to 
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200 400 600 800 1000 

Frequency (Hz) 



Figure 7. Fourier power spectrum for the radial oscillations of 
the centre of a torus (CMR) with a constant angular momentum 
distribution, when perturbed periodically. Three main peaks are 
apparent. From left to right, the first one is the oscillation indcued 
at v per = 200 Hz, the second is the modified epicyclic frequency 
(p— mode) at v\ = k* = 300 Hz, and the third is in a 3:2 relation 
with «*, at frequency U2 = 450 Hz. 




200 400 600 800 1000 

Frequency (Hz) 



Figure 8. Fourier power spectrum for the radial oscillations of 
the centre of a torus (CMR) with a power law distribution of an- 
gular momentum, l{r) oc r a , when perturbed periodically. Sev- 
eral peaks are apparent. The strongest one is at the perturbation 
frequency, u per = 100 Hz. A harmonic at 200 Hz and a weak sub- 
harmonic at 50 Hz are also visible. The fundamental p-mode is at 
v\ = ft* = 250 Hz, and the weaker first overtone at U2 = 360 Hz. 
The latter two are approximately in a 3:2 ratio. 

mimic the effects of strong gravity for the potential pro- 
duced by the central mass. After building initial conditions 
in hydrostatic equilibrium with a prescribed rotation law 
(given through the distribution of specific angular momen- 
tum within the disc) , global impulsive and localized periodic 
perturbations have been applied, and the disc's response has 
been measured through the time variation of certain quan- 
tities, and their corresponding Fourier transforms. An ideal 
gas equation of state has been used throughout. 

As the size of the torus is increased by filling up the 
effective potential well to a greater degree (see Figures 1, 5 
and 6), we find, by applying impulsive kicks, that the fun- 



damental global mode of radial oscillations monotonically 
decreases from the analytic expression for the local radial 
epicyclic frequency at the position of the centre of the torus. 
This effect is present independently of the assumed distri- 
bution of angular momentum (be it a constant or a power 
law) . That is, a slender torus tends to oscillate as a free par- 
ticle would, and a thicker one does so at lower frequencies. 
This effect is simply due to the increased size of the reso- 
nant cavity, and the density stratification within it which 
affects the propagation of sound waves. It was found analyt- 
ically in calculations performed in the Schwarzschild metric 
with height -integr ated tori for inertial-acoustic modes (p- 
modes) recently iRezzolla, Yoshida, Maccarone fc Zanottil 
l2003al:lRezzolla. Yoshida fc Zanottill2003lJ) . and we confirm 
it here for the more general case of a full torus perturbed 
slightly from equilibrium. 

Upon application of a localized and periodic exter- 
nal perturbation, however, the torus responds in a manner 
which does depend on the distribution of angular momen- 
tum. As previously reported in the context of LMXBs (Pa- 
per I), for ^=cst the fundamental radial mode at frequency 
fc* and the first overtone at 1.5k* are excited (although the 
overtone is noticeably weaker than the fundamental). 

For a power law distribution of angular momentum 
(specifically we used £(r) oc r 01 ), the fundamental at k* 
was always present, as was the first overtone at ~ 1.5fv*. 
Furthermore, a harmonic of the perturbation itself, at 2v per , 
and a subharmonic at v per /2 always appeared. We take this 
latter fact to indicate the presence of non-linear mode cou- 
pling, since the applied perturbation was purely sinusoidal 
and contained only a single frequency. 

Trial calculations in both cases revealed that varying 
the perturbation frequency v per by a factor of a few had no 
discernible effect on the results (i.e., this is not the result 
of a coincidence of fundamental frequencies and overtones 
with the particular choice for the perturbation frequency). 

Similar results concerning the relative amplitudes of 
various modes when £(r) is varied can actually be extracted 
from the stu dy of Rezzolla and collaborators. In particular, 
Figure 3 in iRezzolla. Yoshida fc Zanot ti 2003b) shows the 
response (in the density) of a torus with constant angular 
momentum when given an impulsive perturbation. The first 
overtone is weaker than the fundamental by a factor of ~ 30. 
However, results for tori with power-law distributions in £(r) 
show a much larger contrast, with the fundamental and its 
harmonics dominating the frequency spectrum, and p-mode 
overtones being much weaker (about three orders of mag- 
nitude , see Figure 4 in IZanotti. Font. Rezzolla fc Monterq 
(2005)). To complicate matters, it would appear that the 
way in which the pertrubation is applied (e.g., only in the 
density vs. only in the velocity, or both) has an impor- 
ta nt effect on the amplitude of the response (see Figure 5 
in IZanotti. Font. Rezzolla fc Monterol i2005l) ). In calculat- 
ing the eigenfunctions for the height -integrated problem, 
Rezzolla and collaborators have found that these become 
vanishingly small near the inner edge of the torus in certain 
cases, thus making these modes harder to excite. In all cases, 
and as one would expect, they report that the most effi- 
cient way to excite oscillations is to use a perturbation which 
matches the corresponding eigenfunction (i.e., this is when 
the strongest response occurs). We have only attempted one 
type of pertrubation here, given in equation (25), and so 
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have not explored this effect. The vanishing of the eigen- 
function may explain formally why certain oscillations are 
suppressed, or equivalently, why they are so hard to excite 
with arbitrary perturbations. 

In the context of QPOs in LMXBs, producing variabil- 
ity in the X-ray lightcurve from the oscillatio ns described 
here is beyond the scope of this work (see e.g.. ISchnittmanl 
120051 . for how this may actually come about). However, if it 
is actually related to the QPOs, one could argue that the 
presence of definite frequencies, and their relation to one 
another, may serve as an indicator of the distribution of an- 
gular momentum within the disc, or of its relative thickness. 

We believe the study presented here captures the essen- 
tials of the simple system at hand, namely, one with pure 
hydrodynamics in the field of a central object, neglecting self 
gravity and assuming azimuthal symmetry. Clearly in a real 
astrophysical situation these conditions will not be strictly 
met. Magnetic fields are likely to play an important role and 
make the system more complex by introducing additional 
possible modes, as well as coupling to the hydrodynamics. 
Even without magnetic fields, the tori described here will 
not be isolated in a vaccuum, but surrounded by a gaseous 
envelope, the rest of the disk. Mode coupling and/or leakage 
at the interface is likely to affect the oscillatory behaviour 
of the system and is a consideration that deserves further 
study (for a preliminary study addressing these matters see 
lFragildl2005h . 
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